{"cells":[{"cell_type":"markdown","metadata":{"collapsed":false,"id":"185703E2CC3049768049E337F9807E39","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["# Lecture 13: MLM & GLMM\n","\n","## Instructor： 胡传鹏（博士）[Dr. Hu Chuan-Peng]\n","\n","### 南京师范大学心理学院[School of Psychology, Nanjing Normal University]"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"34C89325A5B543E89881E2F4DA2D4EE8","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 回顾：广义线性模型GLM"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"DA5B7512DC6D49F1B2E7AD0AF3C923EB","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["**当因变量为离散变量**\n","\n","比如，因变量为二分变量 (答题正确率)，其中1代表回答正确，0代表回答错误。\n","\n","正确率为离散变量，服从伯努利(Bernoulli)分布。\n","\n","![Image Name](https://cdn.kesci.com/upload/image/rloa62fn8w.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"68E1B4A760654B059BB3408DBEE37861","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["对于因变量为离散变量的情况，我们需要使用广义线性模型(Generalized linear model，GLM)。\n","\n","其特点为：\n","- 分布簇 (dist)不再局限于正态分布，而是允许其他不同的分布，比如 $y \\sim Bernoulli(p)$。\n","- 需要 **链接函数g()** 将 $\\alpha + \\beta * x$  映射到 p 所在的范围。\n","\n","| 一般线性模型 | 广义线性模型 | \n","|---|---|\n","| $y \\sim Normal(\\mu,sigma)$ | $y \\sim dist(p)$ |\n","| $\\mu = \\alpha + \\beta *x$ | $p = g(\\mu)$|"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A1279291B2AD43539815FBD07D0FF622","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["链接函数的具体转化过程，以逻辑(logit)回归为例：\n","1. 令 $z = \\alpha + \\beta *x$，$\\mu$的范围为 $(-\\infty, +\\infty)$。\n","2. $p = g(z)$，其中 g() 为链接函数，输出结果 p 的范围为 $(0,1)$。\n","3.  最后将 p 输入到分布函数中 $y \\sim Bernoulli(p)$。\n","\n","\n","![Image Name](https://cdn.kesci.com/upload/image/rloa6zyf5a.png?imageView2/0/w/600/h/600)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D2A72F5435854179BA7172C4B2DB9AC4","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["## Part 1: HLM&GLMM"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E154452297264C04ACAA90CB14E19868","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["本节的目的在与：了解如何贝叶斯分层广义线性模型 (Bayesian Multilevel Generalized linear model)如何解释数据的分层结构。\n","\n","重点在于：\n","- 了解分层线性模型的基本概念：包括变化截距与变化斜率。\n","- 通过 PyMc 分别拟合变化截距模型与变化斜率模型(包含变化截距)。\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"57C622DB93044A8C93F778338BA1E5F4","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 层级结构\n","#### 我们想用工龄预测某个拥有4个学院的学校员工的工资\n","##### 收集到的数据如下：\n","\n","![Image Name](https://cdn.kesci.com/upload/image/rmcutay7yg.png?imageView2/0/w/960/h/960)\n","\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"0BC5A556E17C487B84965C808E4A49AE","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://cdn.kesci.com/upload/image/rmcuwtjhi2.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D6CBD5DB14FE48D6A3A5C86514FCA8BA","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["将上述例子用广义线性模型来进行表示：\n","$y \\sim Normal(\\mu, \\sigma)$\n","$\\mu  = \\alpha + \\beta *x_{years}$"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"8275A1A13C1241B89353B750FACAECFC","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["由于该数据是一个层级结构，不同的个体可能从属于不同的部门，不同部门间存在的差异可能会对结果的预测产生影响。\n","\n","基于不同的假定，我们可以有不同的对层级进行处理的方式。\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"B4ADC3CBA98849569898CA2E05A81274","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 非池化参数(no pooling)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"80B43C012A994159A23CD4ED8592B15B","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["我们可以假定各部门之间完全没有关联，不同部门之间彼此独立；\n","\n","非池化参数就是各组之间的参数没有关系；\n","\n","在这种情况下，该模型等价于每个部门各有一个广义线性模型。\n","\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"92BEE829B0FF46099FEE7F82B6E8C3C5","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["即：\n","\n","部门1：$y \\sim Normal(\\mu_1, \\sigma)$\n","$\\mu  = \\alpha_1 + \\beta_1 *x_{years}$\n","\n","部门2：$y \\sim Normal(\\mu_1, \\sigma)$\n","$\\mu  = \\alpha_1 + \\beta_1 *x_{years}$\n","\n","...\n","\n","部门j：$y \\sim Normal(\\mu_j, \\sigma)$\n","$\\mu  = \\alpha_j + \\beta_j *x_{years}$\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C02BDF3AA3BA467E99642A9CC255A359","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://cdn.kesci.com/upload/image/rmcwluwquv.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"09F8148B4661451987CDD9CCF457D7D1","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 完全池化参数(complete pooling)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"1D386B05DA2F48A8837960DB77E18BAD","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["我们也可以假定各部门之间完全没有差异；\n","\n","完全池化参数就是各部门间共享同一个参数；\n","\n","在这种情况下，该模型等价于每个部门共享一个广义线性模型。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"896E026B61864399ACCF833A671B9174","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["即：\n","\n","部门1：$y \\sim Normal(\\mu, \\sigma)$\n","$\\mu  = \\alpha + \\beta *x_{years}$\n","\n","部门2：$y \\sim Normal(\\mu, \\sigma)$\n","$\\mu  = \\alpha + \\beta *x_{years}$\n","\n","...\n","\n","部门j：$y \\sim Normal(\\mu, \\sigma)$\n","$\\mu  = \\alpha + \\beta *x_{years}$\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"35875F52A1A44AECB27D800708C3B9F0","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://cdn.kesci.com/upload/image/rmcx1f4dlx.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F1C948B1D6F64664921757E11A5F46E0","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 部分池化参数(partial pooling)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"2FFB5B748A364AC59389BC516B360C78","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["在分层模型里，我们假定各部门之间存在一定的联系\n","\n","部分池化参数就是各组之间并不共享固定的参数值，而是共享用于生成该参数值的同一个概率分布。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"7AD5FC4A39944145BCB2F50366912324","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["即：\n","\n","Group level: $y_{i,j} \\sim Normal(\\mu_{i,j}, \\sigma_{i,j})$\n","$\\mu_{i,j}  = \\alpha_{j} + \\beta_{j} *x_{years}$\n","\n","Population level:$\\alpha_{j} \\sim Normal(\\alpha, \\sigma)$\n","$\\beta_{j} \\sim Normal(\\beta, \\sigma)$"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"DEE1C2681EC3437C85EF717BE5DFECE0","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["在上述形式中，由于线性模型的参数就是它的截距和斜率，所以再更一般的形式中，我们可以直接用$\\theta$表示"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D7CE5ED857CD44549856A8D6F4E07D3B","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://cdn.kesci.com/upload/image/rmcyhbglbb.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6B4B09BA508A4F679B69295DE7673341","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["变化截距与变化斜率会随着不同部门变化，表达式如下：\n","\n","![Image Name](https://cdn.kesci.com/upload/image/rmdh2ogd7u.png?imageView2/0/w/640/h/640)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"8E6A15B6E8C94A24BC00C0D44059964F","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["员工工作年限 (experience)对员工工资 (salary)的预测作用受不同部门的影响\n","\n","- 不同部门员工工资存在不同。\n","\t- 这被称为变化截距，截距的变化代表不同部门间员工平均工资的差异。\n","\t- ![Image Name](https://cdn.kesci.com/upload/image/rmanj9sw62.png?imageView2/0/w/640/h/640)\n","- 不同学校中，员工工作年限对工资的预测作用会发生变化。\n","\t- 这被称为变化斜率，斜率在不同部门间的变化代表，工作年限对工资的预测作用会在不同部门将发生变化。\n","\t- ![Image Name](https://cdn.kesci.com/upload/image/rmank74ytm.png?imageView2/0/w/640/h/640)\n","- 变化截距 + 变化斜率\n","\t- ![Image Name](https://cdn.kesci.com/upload/image/rmankk5jgh.png?imageView2/0/w/640/h/640)\n","\n","图片来源：http://mfviz.com/hierarchical-models/"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"59F165DEFB1B4FB9AB1E6395D3C377B8","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["结合上节课的GLM和这节课的HLM，若因变量为离散变量，数据为嵌套结构时，我们可以讲两者进行结合，这就是GLMM."]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"3236837F73DE4706A82746357EB833B9","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://cdn.kesci.com/upload/image/rmeb2c6lqc.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"code","execution_count":2,"metadata":{"collapsed":false,"id":"337D50C9271F4C1C929F9420B2B593E6","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["WARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n"]}],"source":["#加载需要使用的库\n","%matplotlib inline\n","import numpy as np \n","from scipy import stats\n","import matplotlib.pyplot as plt\n","import pandas as pd\n","import arviz as az\n","import pymc3 as pm\n","\n","np.random.seed(123)  # 变化数种子，确保随后生成的变化数相同"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"DB661350E7094BD8AA327716836B6376","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["## Part 2: 示例1——调查数据"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"86222DC0EDEC4EE2BDA4B0FA9BEE42C3","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["**具体实例**\n","\n","基于泰国初等教育数据库，我们想探究，学生是否接受学前教育对之后学生是否留级的影响？"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9FC53DA0FF12486095EBF73F4D1E4C11","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["数据来源：泰国初等教育的全国调查 (Raudenbush & Bhumirat，1992)。\n","- 数据中的每一行表示一名学生。\n","- 因变量 `REPEAT` 为二分类变量，表示学生在初等教育期间是否留级，1 = 留级，0 = 非留级。\n","- 自变量为学生是否受过学前教育 `PPED` (0 = 否，1 = 是)。\n","\n","数据来源：https://github.com/MultiLevelAnalysis/Datasets-third-edition-Multilevel-book/blob/master/chapter%206/Thaieduc/thaieduc.sav"]},{"cell_type":"code","execution_count":3,"metadata":{"collapsed":false,"id":"CA15DD19228048AF89E9E4785D94425C","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["# 加载数据\n","data = pd.read_csv(\"/home/mw/input/data1317/thaieduc.csv\")\n","data.PPED = data.PPED.map({'yes':1,'no':0})"]},{"cell_type":"code","execution_count":4,"metadata":{"collapsed":false,"id":"A50A3E682CEB4E95892F664D4DACD83E","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>SCHOOLID</th>\n","      <th>SEX</th>\n","      <th>PPED</th>\n","      <th>REPEAT</th>\n","      <th>MSESC</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>0</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>0</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["   SCHOOLID   SEX  PPED  REPEAT  MSESC\n","0     10103  girl     0       0   0.88\n","1     10103  girl     0       0   0.88\n","2     10103  girl     1       0   0.88\n","3     10103  girl     1       0   0.88\n","4     10103  girl     1       0   0.88"]},"execution_count":4,"metadata":{},"output_type":"execute_result"}],"source":["data.head()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"8CEA3D451F0E432BA2B2CD63E4B6F1BB","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### Workflow\n","\n","在回顾广义线性模型 (GLM)与介绍分层线性模型 (MLM/HLM)的基础知识后，我们回到泰国初等教育的实例，并通过PyMC完成的贝叶斯分层广义线性模型建模的全过程 (full workflow)。\n","\n","![Image Name](https://cdn.kesci.com/upload/image/rkvikqg9q6.png?imageView2/0/w/650/h/650)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9058B214F6F6413AB1F6CAA9770D2B0C","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (1) 提出研究问题\n","\n","使用泰国初等教育数据探究的研究问题为：\n","- 学生是否接受学前教育对学生留级的影响？\n","- 在考虑不同学校的分层结构时，接受学前教育对学生是否留级的影响有哪些变化？\n","- 分层模型是否能提供额外的信息？"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A0057168DAEA4752A27834F00113CB29","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (2) 数据收集"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"00744E1AC20248219ACBD22038ADBBE6","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["数据来源：泰国初等教育的全国调查 (Raudenbush & Bhumirat，1992)。\n","- 数据中的每一行表示一名学生。\n","- 因变量 `REPEAT` 为二分类变量，表示学生在初等教育期间是否留级，1 = 留级，0 = 非留级。\n","- 自变量为学生是否受过学前教育 `PPED` (0 = 否，1 = 是)。\n","- 其他变量包括学生性别 `SEX`、学生所在的学校 `SCHOOLID`和学校平均社会经济地位 (SES)分数 `MSESC`。\n","\n","数据来源：https://github.com/MultiLevelAnalysis/Datasets-third-edition-Multilevel-book/blob/master/chapter%206/Thaieduc/thaieduc.sav"]},{"cell_type":"code","execution_count":5,"metadata":{"collapsed":false,"id":"EB626C0FB2B74C50B3A5A8FDEA6F81DF","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["# 加载数据\n","data = pd.read_csv(\"/home/mw/input/data1317/thaieduc.csv\")\n","data.PPED = data.PPED.map({'yes':1,'no':0})"]},{"cell_type":"code","execution_count":6,"metadata":{"collapsed":false,"id":"750641CE49A24B8CA2047E8D5FD62F94","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>SCHOOLID</th>\n","      <th>SEX</th>\n","      <th>PPED</th>\n","      <th>REPEAT</th>\n","      <th>MSESC</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>0</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>0</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["   SCHOOLID   SEX  PPED  REPEAT  MSESC\n","0     10103  girl     0       0   0.88\n","1     10103  girl     0       0   0.88\n","2     10103  girl     1       0   0.88\n","3     10103  girl     1       0   0.88\n","4     10103  girl     1       0   0.88"]},"execution_count":6,"metadata":{},"output_type":"execute_result"}],"source":["data.head()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"B41C5FF34A5E46C0A49B0FABA2A6BD35","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 描述统计与探索性分析"]},{"cell_type":"code","execution_count":7,"metadata":{"collapsed":false,"id":"61ED9314257F4A99B162E1805219BFF8","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["PPED\n","0    0.277778\n","1    0.191729\n","Name: REPEAT, dtype: float64"]},"execution_count":7,"metadata":{},"output_type":"execute_result"}],"source":["# groupby 可以按照某条件对数据进行分组，之后再使用 mean函数可以实现分组求均值。\n","data.groupby('PPED').REPEAT.mean() "]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"1AEFD54BF1554DC89F6F6258DB700422","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["可以发现，没有接受过学前教育学生的留级率 (18%)高于接受过学前教育的学生 (11%)。\n","\n","可视化如下："]},{"cell_type":"code","execution_count":8,"metadata":{"collapsed":false,"id":"AF46B943652441D99746926AD21849E7","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/F00500D64509449992C588D3DFE24B1C/rmd0j7p0u3.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["data.groupby('PPED').REPEAT.mean().plot.bar()\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E2C3ADABBC3F40A38C81CB27492DF4AA","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["对比不同学校下，学前教育对留级率的影响。\n","\n","其中，每一条线代表一个学校。\n","\n","可见，不同学校间的学前教育对留级率的影响存在差异。但整体存在接受学前教育后留级率下降的趋势。"]},{"cell_type":"code","execution_count":9,"metadata":{"collapsed":false,"id":"56D50697EC5B4211B562D67644364489","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["Text(0, 0.5, 'REPEAT \\\\%')"]},"execution_count":9,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/B1BE3856CCAD443D82F39C0239E3B8A2/rmd0jm3rav.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 计算不同学校下，是否接受学前教育的留级率\n","data_plot = data.groupby([\"SCHOOLID\",'PPED']).REPEAT.mean() \n","# 根据留级率进行绘图\n","for schoolID,data_i in data_plot.groupby(\"SCHOOLID\"):\n","    data_i.plot.line(x=\"PPED\")\n","plt.xlabel(\"pped\")             # 横坐标为是否接受学前教育\n","plt.xticks([0,1],[\"no\",\"yes\"]) # 0为否，1为是\n","plt.ylabel(r\"REPEAT \\%\")       # 纵坐标为留级率"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"950D532F6E7243CDAFFDF998C43BC02F","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (3) 选择模型"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D4E3C26734734308B98CD5562D5DAD22","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["由于因变量(是否留级)是二分离散变量，因此我们选择使用基于伯努利(Bernoulli)分布的广义线性模型(Generalized linear model，GLM)进行模型拟合。\n","\n","这里我们考虑三种可能的模型：\n","- 模型1：仅包含变化截距的模型。表示拟合不同学校的平均留级率。\n","- 模型2：变化截距模型。在模型1的基础上，考虑学前教育在level1 (group)的总体效应，这种总体效应也成为固定效应。\n","- 模型3：变化斜率模型。在模型2的基础上，考虑学前教育在level2 (population)的不同学校上的不同效应，这种总体效应也成为变化效应。变化截距也属于变化效应。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"023ECC7CEBFE42898BCEB6693B68EE46","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 模型1：仅包含变化截距的模型"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"CF464C99DC1F4A1FAF9177ADCAEA7351","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["$$\n","y_{(i,j)} \\sim Bernoulli(p_{(i,j)})\n","$$\n","该模型只拟合不同学校的平均留级率。\n","- 其中i为不同学生，j为不同学校。\n","\n","$$\n","p_{(i,j)} = g(\\alpha_j) \\\\\n","\\alpha_j \\sim Normal(\\mu_{\\alpha},\\sigma_{\\alpha})\n","$$\n","- $\\alpha_{j}$ 为截距，代表平均留级率。该解决随着不同学校j进行变化，因此称为变化截距。\n","- g()为链接函数\n","- 其中，$\\mu_{\\alpha}, \\sigma_{\\alpha}$ 为 level 1 group 的效应，而 $\\alpha_{j}$ 为 level 2 population 的效应。\n","- 注意，此时的模型没有考虑自变量 **学前教育** 的影响。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"BFD8636763B44F4482B43D3E0253E945","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["**先验**\n","\n","此时模型参数 ($\\alpha$)的先验与 $\\mu_{\\alpha}, \\sigma_{\\alpha}$有关。而 $\\mu_{\\alpha}, \\sigma_{\\alpha}$ 也服从于不同的分布。比如：\n","\n","$$\n","\\mu_{\\alpha} \\sim Normal(0,100) \\\\ \n","\\sigma_{\\alpha} \\sim HalfNormal(50)\n","$$\n","- 其中 $\\mu_{\\alpha}, \\sigma_{\\alpha}$ 称为超参 (hyperparameters)，可以理解为超越先验 $\\alpha$ 的参数。\n","- $Normal(0,100),HalfNormal(50)$ 称为超先验 (hyperprior)，即先验的先验。"]},{"cell_type":"code","execution_count":10,"metadata":{"collapsed":false,"id":"A3D0827E19084130AAD0CDBE2FE451A6","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["# 将数据分层变量\"学校(school)\"转换为因子(factor)类型\n","school_idxs, school = pd.factorize(data.SCHOOLID)\n","# 定义学校与数据的映射：即标注哪名学生(行)属于哪一所学校\n","coords = {\n","    \"school\": school,\n","    \"obs_id\": np.arange(len(school_idxs)),\n","}"]},{"cell_type":"code","execution_count":11,"metadata":{"collapsed":false,"id":"651223ECD75D4E3D9DA91A1F1EA1D32F","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model(coords=coords) as model1:\n","    # 定义超先验 Hyperpriors\n","    mu_alpha = pm.Normal(\"mu_alpha\", mu=0.0, sigma=100)  # 对应上述公式\n","    sigma_alpha = pm.HalfNormal(\"sigma_alpha\", 50)\n","    # 定义先验\n","    alpha = pm.Normal('alpha',mu=mu_alpha,sd=sigma_alpha, dims=\"school\")\n","    # 定义数据分层变量\"学校(school)\"\n","    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n","    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n","    p = pm.Deterministic(\"p\", alpha[school_idx])\n","    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"399F11566DFF4AE6B269A485F9B5C801","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["我们通过 PyMc 自带的函数 `model_to_graphviz()` 将模型可视化。\n","\n","需要注意的是：\n","- 我们假设 $\\alpha$ (先验)服从一个正态分布 其均值为 mu_alpha，误差为 sigma_alpha。\n","- 并且 mu_alpha 和 sigma_alpha 各自服从不同的分布，如下图，这些分布为 **超先验(hyperpriors)**。\n","- 定义超先验的目的在于约束学校变量带来的差异，这被称为 \"shrinkage\"，是分层模型 partial pooling的效果。更多详情请参考 https://mc-stan.org/rstanarm/articles/pooling.html"]},{"cell_type":"code","execution_count":12,"metadata":{"collapsed":false,"id":"A44367B7A4794E2BAF94D5D4808D694C","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/7B7DF272235145FE85164CFA7EF4340D/rmd0ll1c10.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7fe21067dd10>"]},"execution_count":12,"metadata":{},"output_type":"execute_result"}],"source":["pm.model_to_graphviz(model1)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E625BE6E4B5D46698820864A760EE460","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 模型2：变化截距模型\n","\n","在模型1的基础上，考虑学前教育在level1的总体效应。\n","$$\n","y_{(i,j)} \\sim Bernoulli(p_{(i,j)}) \\\\\n","p_{(i,j)} = g(\\alpha_j + \\beta*\\text{x}_{i, j}) \\\\\n","\\alpha_j \\sim Normal(\\mu_{\\alpha},\\sigma_{\\alpha})\n","$$\n","- x为自变量是否接受学前教育。\n","- 需要注意的是，学前教育的效应 $\\beta$ 不随学校进行变化。这意味着它在不同学生和不同学校间保持固定，因此也称为**固定效应**。\n","\n","\n","**先验**\n","\n","$$\n","\\mu_{\\alpha} \\sim Normal(0,100) \\\\ \n","\\sigma_{\\alpha} \\sim HalfNormal(50) \\\\\n","\\beta \\sim Normal(0,100) \\\\ \n","$$"]},{"cell_type":"code","execution_count":13,"metadata":{"collapsed":false,"id":"5712F817B57F4C46919ACA11D11ADD43","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model(coords=coords) as model2:\n","    # 定义超先验 Hyperpriors\n","    mu_alpha = pm.Normal(\"mu_alpha\", mu=0.0, sigma=100)\n","    sigma_alpha = pm.HalfNormal(\"sigma_alpha\", 50)\n","    # 定义先验\n","    alpha = pm.Normal('alpha',mu=mu_alpha,sd=sigma_alpha, dims=\"school\")\n","    # x为自变量，是之前已经载入的数据\n","    x = pm.Data(\"x\", data.PPED, dims=\"obs_id\")\n","    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n","   \n","    #####################################\n","    # 与模型1的区别在于定义固定效应 beta, 以及线性公式p的变化\n","    #####################################\n","    beta = pm.Normal('beta',mu=0,sd=10)\n","    p = pm.Deterministic(\"p\", alpha[school_idx] + beta * x)\n","    \n","    # Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")"]},{"cell_type":"code","execution_count":14,"metadata":{"collapsed":false,"id":"BF58C55B0156463D8A5EB7CA7B52BE70","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/A11D835D5928487885D6EE4993ACA506/rmd0m5cjch.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7fe1fea42bd0>"]},"execution_count":14,"metadata":{},"output_type":"execute_result"}],"source":["pm.model_to_graphviz(model2)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E420DE3E91574187AC94EDDE8AE87C96","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 模型3：变化截距+斜率模型\n","\n","$$\n","y_{(i,j)} \\sim Bernoulli(p_{(i,j)})\n","p_{(i,j)} = g(\\alpha_j + \\beta*\\text{x}_{i, j}) \\\\\n","\\alpha_j \\sim Normal(\\mu_{\\alpha},\\sigma_{\\alpha}) \\\\\n","\\beta_j \\sim Normal(\\mu_{\\beta},\\sigma_{\\beta})\n","$$\n","在模型2的基础上，考虑学前教育在level2的不同学校上的不同效\n","- 此时，x的效应 $\\beta_{j}$ 随不同的学校 j 进行变化。\n","- 这种总体效应也成为变化效应。其实变化截距也属于变化效应。\n","\n","\n","**先验**\n","\n","$$\n","\\mu_{\\alpha} \\sim Normal(0,100) \\\\ \n","\\sigma_{\\alpha} \\sim HalfNormal(50) \\\\\n","\\mu_{\\beta} \\sim Normal(0,100) \\\\ \n","\\sigma_{\\beta} \\sim HalfNormal(50) \\\\\n","$$"]},{"cell_type":"code","execution_count":16,"metadata":{"collapsed":false,"id":"87BE6D7693AB40B98483132190588582","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model(coords=coords) as model3:\n","    #####################################\n","    # 与模型2的区别在于，额外定义效应beta的超先验 Hyperpriors\n","    #####################################\n","    mu_alpha = pm.Normal(\"mu_alpha\", mu=0.0, sigma=100)\n","    sigma_alpha = pm.HalfNormal(\"sigma_alpha\", 50)\n","    mu_beta = pm.Normal(\"mu_beta\", mu=0.0, sigma=100)\n","    sigma_beta = pm.HalfNormal(\"sigma_beta\", 50)\n","    \n","    # 定义先验\n","    alpha = pm.Normal('alpha',mu=mu_alpha,sd=sigma_alpha, dims=\"school\")\n","    beta = pm.Normal('beta',mu=mu_beta,sd=sigma_beta, dims=\"school\")\n","    # x为自变量，是之前已经载入的数据\n","    x = pm.Data(\"x\", data.PPED, dims=\"obs_id\")\n","    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n","    \n","    #####################################\n","    # 与模型2的区别在于，注意此时的beta随着学校school_idx进行变化\n","    #####################################\n","    p = pm.Deterministic(\"p\", alpha[school_idx] + beta[school_idx] * x)\n","    # Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"90D3D7BE9535483B9B829E2DD408CD49","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["注意观察，每个变化效应的先验 (alpha 和 beta)都存在超先验 (hyperpriors)。\n","\n","同样，beta参数的超先验约束了不同学校间学前教育效应的影响，可避免出现极端值。"]},{"cell_type":"code","execution_count":17,"metadata":{"collapsed":false,"id":"B18204031F2E46C89339BE4E30675A08","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/7EC31CF87DEC49578B8884E0D621BE66/rmd0na2u4p.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7fe1f632dd50>"]},"execution_count":17,"metadata":{},"output_type":"execute_result"}],"source":["pm.model_to_graphviz(model3)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"B8F626D93A9349699E270D12294F980A","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["这里我们先考虑最后一个模型 (模型3)，在后面模型比较的部分，再考虑其他两个模型。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"8C1DFFE75A1C43F9B741EA14B2B1DE7C","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (4)选择先验"]},{"cell_type":"code","execution_count":18,"metadata":{"collapsed":false,"id":"D0D419B5A51D453A9FA1D3F1A6CAB89D","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with model3:\n","    # 先验预测检查\n","    prior_checks = pm.sample_prior_predictive(samples=500)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"8FD35F2C8E7442049529C89E90685137","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["首先查看 p 的先验分布。\n","\n","注意，这里没有用链接函数对 p 进行转换，因此其范围在 -200到200之间，而不是 0-1之间。"]},{"cell_type":"code","execution_count":19,"metadata":{"collapsed":false,"id":"8C53E250569544B5B55E16AAF4D11097","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/data/base.py:174: UserWarning: More chains (500) than draws (338). Passed array should have shape (chains, draws, *shape)\n","  UserWarning,\n"]},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/17AF2F350747449290D1901B6B8DF781/rmd0of8zl6.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_density(\n","    {'p':prior_checks['p']}\n","    )\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"26BD0368929B42A08505CFCA059BA6ED","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["其次，查看模型截距和斜率的先验。"]},{"cell_type":"code","execution_count":20,"metadata":{"collapsed":false,"id":"A169D5188F424783B18565E00325F3B1","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/data/base.py:174: UserWarning: More chains (500) than draws (19). Passed array should have shape (chains, draws, *shape)\n","  UserWarning,\n"]},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/135F1944152A43EF80602D89B46A3512/rmd0orxz7l.png\">"],"text/plain":["<Figure size 993.6x331.2 with 2 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_density(\n","    {'alpha':prior_checks['alpha'],\n","    'beta':prior_checks['beta']}\n","    )\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D86950D7192C4B3CA6E3CFF3A04C4B5F","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["最后，别忘了查看超先验所对应超参数(hyperparameters)的先验分布。"]},{"cell_type":"code","execution_count":21,"metadata":{"collapsed":false,"id":"A9B8A6C9ADC1461CA8262E498B43DAED","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/0DA60DC5A40D4A53A4C6641777F2A68F/rmd0p6g3ea.png\">"],"text/plain":["<Figure size 1490.4x662.4 with 6 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_density(\n","    {\n","        'mu_alpha':prior_checks['mu_alpha'],\n","        'sigma_alpha':prior_checks['sigma_alpha'],\n","        'mu_beta':prior_checks['mu_beta'],\n","        'sigma_beta':prior_checks['sigma_beta']\n","    }\n","    )\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F410BA1D9283474B859DE5333FD03E1E","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (5) 拟合数据"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F70FB11DC45C4032A2656206272B9963","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["拟合数据需要注意，虽然我们主要探究的自变量为\"是否接受学前教育\"，但是数据分层结构变量\"学校\"也是非常重要的，不要忽视了。"]},{"cell_type":"code","execution_count":22,"metadata":{"collapsed":false,"id":"6B3715000B214782A29DF3C8996F8974","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>SCHOOLID</th>\n","      <th>SEX</th>\n","      <th>PPED</th>\n","      <th>REPEAT</th>\n","      <th>MSESC</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>0</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>0</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>10103</td>\n","      <td>girl</td>\n","      <td>1</td>\n","      <td>0</td>\n","      <td>0.88</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["   SCHOOLID   SEX  PPED  REPEAT  MSESC\n","0     10103  girl     0       0   0.88\n","1     10103  girl     0       0   0.88\n","2     10103  girl     1       0   0.88\n","3     10103  girl     1       0   0.88\n","4     10103  girl     1       0   0.88"]},"execution_count":22,"metadata":{},"output_type":"execute_result"}],"source":["data.head()"]},{"cell_type":"code","execution_count":23,"metadata":{"collapsed":false,"id":"D0A5D584A4BC4E23B8AC130990049623","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/B519FED7525546C88CD6974D1AA5246F/rmd0pzkoq5.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["data.groupby('SCHOOLID').REPEAT.mean().plot.bar()\n","plt.show()"]},{"cell_type":"code","execution_count":24,"metadata":{"collapsed":false,"id":"7EC1A4B9F2794777963C149ECBB5CC48","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["# 将数据分层变量\"学校(school)\"转换为因子(factor)类型\n","school_idxs, school = pd.factorize(data.SCHOOLID)\n","# 定义学校与数据的映射：即标注哪名学生(行)属于哪一所学校\n","coords = {\n","    \"school\": school,\n","    \"obs_id\": np.arange(len(school_idxs)),\n","}\n","\n","with pm.Model(coords=coords) as model3:\n","    # 定义level2学校的 Hyperpriors\n","    mu_alpha = pm.Normal(\"mu_alpha\", mu=0.0, sigma=100)\n","    sigma_alpha = pm.HalfNormal(\"sigma_alpha\", 50)\n","    mu_beta = pm.Normal(\"mu_beta\", mu=0.0, sigma=100)\n","    sigma_beta = pm.HalfNormal(\"sigma_beta\", 50)\n","    # 定义先验\n","    alpha = pm.Normal('alpha',mu=mu_alpha,sd=sigma_alpha, dims=\"school\")\n","    beta = pm.Normal('beta',mu=mu_beta,sd=sigma_beta, dims=\"school\")\n","    # x为自变量，是之前已经载入的数据\n","    x = pm.Data(\"x\", data.PPED, dims=\"obs_id\")\n","    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n","    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n","    p = pm.Deterministic(\"p\", alpha[school_idx] + beta[school_idx] * x)\n","    # Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6E899D1D69314DB8BA8BC73CFA6303D2","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (6)采样过程诊断\n","\n","如果使用MCMC对后验进行近似，则需要首先对MCMC过程进行评估。\n","\n","* 是否收敛；\n","* 是否接近真实的后验。\n","\n","对采样过程的评估我们会采用目视检查或rhat这个指标"]},{"cell_type":"code","execution_count":25,"metadata":{"collapsed":false,"id":"71487D6AB3A04A0DBEB550ADEEA71B3C","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [beta, alpha, sigma_beta, mu_beta, sigma_alpha, mu_alpha]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 22 seconds.\n"]},{"name":"stdout","output_type":"stream","text":["0, dim: obs_id, 338 =? 338\n"]},{"name":"stderr","output_type":"stream","text":["There were 29 divergences after tuning. Increase `target_accept` or reparameterize.\n","There were 339 divergences after tuning. Increase `target_accept` or reparameterize.\n","The acceptance probability does not match the target. It is 0.752539189964468, but should be close to 0.9. Try to increase the number of tuning steps.\n","The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.\n","The estimated number of effective samples is smaller than 200 for some parameters.\n"]}],"source":["with model3:\n","    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n","    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n","    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n","    trace3 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"code","execution_count":26,"metadata":{"collapsed":false,"id":"8201E58CAEF142EA80F22A5DC3885EE2","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/4C48AC88E0A3460F8D4C5435B82EDBDB/rmd0spg432.png\">"],"text/plain":["<Figure size 864x288 with 4 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_trace(trace3, var_names=['alpha','beta'])\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C093677234774F02A3B7B142DB5D74D1","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["注意，由于参数会随着学校进行变化，因此每一所学校(19所学校)都对应一个参数，如下："]},{"cell_type":"code","execution_count":27,"metadata":{"collapsed":false,"id":"CA00B4F785E44B788D758B309EB5D0DF","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>alpha[0]</th>\n","      <td>0.029</td>\n","      <td>0.021</td>\n","      <td>323.0</td>\n","      <td>1227.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[1]</th>\n","      <td>0.014</td>\n","      <td>0.010</td>\n","      <td>850.0</td>\n","      <td>1607.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[2]</th>\n","      <td>0.019</td>\n","      <td>0.013</td>\n","      <td>611.0</td>\n","      <td>1257.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[3]</th>\n","      <td>0.025</td>\n","      <td>0.018</td>\n","      <td>447.0</td>\n","      <td>1065.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[4]</th>\n","      <td>0.035</td>\n","      <td>0.025</td>\n","      <td>251.0</td>\n","      <td>1148.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[5]</th>\n","      <td>0.010</td>\n","      <td>0.007</td>\n","      <td>1211.0</td>\n","      <td>1739.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[6]</th>\n","      <td>0.013</td>\n","      <td>0.009</td>\n","      <td>864.0</td>\n","      <td>1423.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[7]</th>\n","      <td>0.035</td>\n","      <td>0.025</td>\n","      <td>261.0</td>\n","      <td>1012.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[8]</th>\n","      <td>0.015</td>\n","      <td>0.011</td>\n","      <td>660.0</td>\n","      <td>321.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[9]</th>\n","      <td>0.017</td>\n","      <td>0.013</td>\n","      <td>731.0</td>\n","      <td>1455.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[10]</th>\n","      <td>0.023</td>\n","      <td>0.017</td>\n","      <td>472.0</td>\n","      <td>1299.0</td>\n","      <td>1.00</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[11]</th>\n","      <td>0.014</td>\n","      <td>0.010</td>\n","      <td>909.0</td>\n","      <td>1382.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[12]</th>\n","      <td>0.022</td>\n","      <td>0.016</td>\n","      <td>523.0</td>\n","      <td>947.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[13]</th>\n","      <td>0.018</td>\n","      <td>0.013</td>\n","      <td>598.0</td>\n","      <td>1009.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[14]</th>\n","      <td>0.023</td>\n","      <td>0.016</td>\n","      <td>494.0</td>\n","      <td>1027.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[15]</th>\n","      <td>0.017</td>\n","      <td>0.012</td>\n","      <td>837.0</td>\n","      <td>1052.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[16]</th>\n","      <td>0.018</td>\n","      <td>0.013</td>\n","      <td>648.0</td>\n","      <td>686.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[17]</th>\n","      <td>0.019</td>\n","      <td>0.015</td>\n","      <td>682.0</td>\n","      <td>1366.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>alpha[18]</th>\n","      <td>0.018</td>\n","      <td>0.013</td>\n","      <td>672.0</td>\n","      <td>1231.0</td>\n","      <td>1.00</td>\n","    </tr>\n","    <tr>\n","      <th>beta[0]</th>\n","      <td>0.048</td>\n","      <td>0.034</td>\n","      <td>199.0</td>\n","      <td>932.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[1]</th>\n","      <td>0.020</td>\n","      <td>0.015</td>\n","      <td>763.0</td>\n","      <td>1690.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[2]</th>\n","      <td>0.019</td>\n","      <td>0.013</td>\n","      <td>842.0</td>\n","      <td>1817.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>beta[3]</th>\n","      <td>0.028</td>\n","      <td>0.020</td>\n","      <td>458.0</td>\n","      <td>1097.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[4]</th>\n","      <td>0.080</td>\n","      <td>0.057</td>\n","      <td>124.0</td>\n","      <td>1198.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[5]</th>\n","      <td>0.094</td>\n","      <td>0.067</td>\n","      <td>100.0</td>\n","      <td>271.0</td>\n","      <td>1.03</td>\n","    </tr>\n","    <tr>\n","      <th>beta[6]</th>\n","      <td>0.040</td>\n","      <td>0.028</td>\n","      <td>284.0</td>\n","      <td>1002.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[7]</th>\n","      <td>0.034</td>\n","      <td>0.024</td>\n","      <td>393.0</td>\n","      <td>1226.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[8]</th>\n","      <td>0.034</td>\n","      <td>0.024</td>\n","      <td>336.0</td>\n","      <td>879.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[9]</th>\n","      <td>0.053</td>\n","      <td>0.038</td>\n","      <td>178.0</td>\n","      <td>1251.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[10]</th>\n","      <td>0.028</td>\n","      <td>0.020</td>\n","      <td>570.0</td>\n","      <td>1433.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[11]</th>\n","      <td>0.053</td>\n","      <td>0.038</td>\n","      <td>138.0</td>\n","      <td>1085.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[12]</th>\n","      <td>0.022</td>\n","      <td>0.016</td>\n","      <td>686.0</td>\n","      <td>1179.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[13]</th>\n","      <td>0.023</td>\n","      <td>0.016</td>\n","      <td>487.0</td>\n","      <td>1561.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[14]</th>\n","      <td>0.047</td>\n","      <td>0.033</td>\n","      <td>160.0</td>\n","      <td>1128.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>beta[15]</th>\n","      <td>0.021</td>\n","      <td>0.016</td>\n","      <td>1121.0</td>\n","      <td>1754.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[16]</th>\n","      <td>0.018</td>\n","      <td>0.013</td>\n","      <td>937.0</td>\n","      <td>1664.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[17]</th>\n","      <td>0.036</td>\n","      <td>0.026</td>\n","      <td>297.0</td>\n","      <td>1435.0</td>\n","      <td>1.00</td>\n","    </tr>\n","    <tr>\n","      <th>beta[18]</th>\n","      <td>0.041</td>\n","      <td>0.029</td>\n","      <td>204.0</td>\n","      <td>1436.0</td>\n","      <td>1.01</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["           mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat\n","alpha[0]       0.029    0.021     323.0    1227.0   1.01\n","alpha[1]       0.014    0.010     850.0    1607.0   1.01\n","alpha[2]       0.019    0.013     611.0    1257.0   1.02\n","alpha[3]       0.025    0.018     447.0    1065.0   1.02\n","alpha[4]       0.035    0.025     251.0    1148.0   1.01\n","alpha[5]       0.010    0.007    1211.0    1739.0   1.01\n","alpha[6]       0.013    0.009     864.0    1423.0   1.02\n","alpha[7]       0.035    0.025     261.0    1012.0   1.01\n","alpha[8]       0.015    0.011     660.0     321.0   1.02\n","alpha[9]       0.017    0.013     731.0    1455.0   1.01\n","alpha[10]      0.023    0.017     472.0    1299.0   1.00\n","alpha[11]      0.014    0.010     909.0    1382.0   1.01\n","alpha[12]      0.022    0.016     523.0     947.0   1.01\n","alpha[13]      0.018    0.013     598.0    1009.0   1.02\n","alpha[14]      0.023    0.016     494.0    1027.0   1.02\n","alpha[15]      0.017    0.012     837.0    1052.0   1.02\n","alpha[16]      0.018    0.013     648.0     686.0   1.01\n","alpha[17]      0.019    0.015     682.0    1366.0   1.01\n","alpha[18]      0.018    0.013     672.0    1231.0   1.00\n","beta[0]        0.048    0.034     199.0     932.0   1.01\n","beta[1]        0.020    0.015     763.0    1690.0   1.01\n","beta[2]        0.019    0.013     842.0    1817.0   1.02\n","beta[3]        0.028    0.020     458.0    1097.0   1.01\n","beta[4]        0.080    0.057     124.0    1198.0   1.01\n","beta[5]        0.094    0.067     100.0     271.0   1.03\n","beta[6]        0.040    0.028     284.0    1002.0   1.01\n","beta[7]        0.034    0.024     393.0    1226.0   1.01\n","beta[8]        0.034    0.024     336.0     879.0   1.01\n","beta[9]        0.053    0.038     178.0    1251.0   1.01\n","beta[10]       0.028    0.020     570.0    1433.0   1.01\n","beta[11]       0.053    0.038     138.0    1085.0   1.01\n","beta[12]       0.022    0.016     686.0    1179.0   1.01\n","beta[13]       0.023    0.016     487.0    1561.0   1.01\n","beta[14]       0.047    0.033     160.0    1128.0   1.02\n","beta[15]       0.021    0.016    1121.0    1754.0   1.01\n","beta[16]       0.018    0.013     937.0    1664.0   1.01\n","beta[17]       0.036    0.026     297.0    1435.0   1.00\n","beta[18]       0.041    0.029     204.0    1436.0   1.01"]},"execution_count":27,"metadata":{},"output_type":"execute_result"}],"source":["az.summary(trace3, var_names=['alpha','beta'], kind=\"diagnostics\")"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C1C6D51DC5C441DE97CEE1C2AD4CE693","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (7)模型诊断\n","\n","在MCMC有效的前提下，需要继续检验模型是否能够较好地拟合数据。\n","\n","我们会使用后验预测分布通过我们得到的参数生成一批模拟数据，并将其与真实数据进行对比。"]},{"cell_type":"code","execution_count":28,"metadata":{"collapsed":false,"id":"7DC2F15DC83E448A95DC976F01C735A7","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["\n","    <div>\n","      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n","      100.00% [4000/4000 00:03&lt;00:00]\n","    </div>\n","    "],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"update_display_data"},{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n","  FutureWarning,\n"]}],"source":["# 后验预测分布的计算仍在容器中进行\n","with model3:\n","    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n","    ppc_y = pm.sample_posterior_predictive(trace3.posterior) \n","#将ppc_y转化为InferenceData对象合并到trace中\n","az.concat(trace3, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)"]},{"cell_type":"code","execution_count":29,"metadata":{"collapsed":false,"id":"C7FF6B69932449F582E01603FCED6273","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n","  fig.canvas.print_figure(bytes_io, **kw)\n"]},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/E1A54B5D7C01421F97C4537CF0AE8878/rmd0t5wtic.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制后验预测分布\n","az.plot_ppc(trace3)\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"68E80371BF9D4358AA5D212F358AC2CB","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (8)模型比较"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E9DBFC0AF2634F61A1B0CBCDAB8B89B5","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["前面的诊断过程，我们只考虑了模型3 (变化截距+斜率模型)。在模型比较阶段，我们可以同时比较三个模型：\n","- 模型1：仅包含变化截距的模型，拟合不同学校的平均留级率。\n","- 模型2：变化截距模型。在模型1的基础上，考虑学前教育在group level的总体效应。\n","- 模型3：变化截距+斜率模型。在模型2的基础上，考虑学前教育在不同学校上 (population level)的不同效应。"]},{"cell_type":"code","execution_count":30,"metadata":{"collapsed":false,"id":"CBFABC8D84EB417A8A6F8BF0B7510353","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [alpha, sigma_alpha, mu_alpha]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 9 seconds.\n","There were 8 divergences after tuning. Increase `target_accept` or reparameterize.\n","The number of effective samples is smaller than 25% for some parameters.\n"]},{"name":"stdout","output_type":"stream","text":["0, dim: obs_id, 338 =? 338\n"]},{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [beta, alpha, sigma_alpha, mu_alpha]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 11 seconds.\n","There were 49 divergences after tuning. Increase `target_accept` or reparameterize.\n","There were 15 divergences after tuning. Increase `target_accept` or reparameterize.\n","The estimated number of effective samples is smaller than 200 for some parameters.\n"]},{"name":"stdout","output_type":"stream","text":["0, dim: obs_id, 338 =? 338\n"]}],"source":["# 对模型进行mcmc采样\n","with model1:\n","    trace1 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)\n","with model2:\n","    trace2 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)\n","# 模型3前面完成了采样，这里可以不用再采样了\n","# with model3:\n","#     trace2 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"code","execution_count":31,"metadata":{"collapsed":false,"id":"C3433A59113149FEB89E1E7570CF7E3A","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n","  \"The default method used to estimate the weights for each model,\"\n"]},{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>rank</th>\n","      <th>loo</th>\n","      <th>p_loo</th>\n","      <th>d_loo</th>\n","      <th>weight</th>\n","      <th>se</th>\n","      <th>dse</th>\n","      <th>warning</th>\n","      <th>loo_scale</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>随机斜率模型</th>\n","      <td>0</td>\n","      <td>-171.263208</td>\n","      <td>12.461399</td>\n","      <td>0.000000</td>\n","      <td>6.657864e-01</td>\n","      <td>10.313712</td>\n","      <td>0.000000</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","    <tr>\n","      <th>仅包含截距的模型</th>\n","      <td>1</td>\n","      <td>-171.747158</td>\n","      <td>10.166437</td>\n","      <td>0.483949</td>\n","      <td>3.342136e-01</td>\n","      <td>10.156738</td>\n","      <td>1.729660</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","    <tr>\n","      <th>随机截距模型</th>\n","      <td>2</td>\n","      <td>-171.944453</td>\n","      <td>10.723415</td>\n","      <td>0.681245</td>\n","      <td>4.440892e-16</td>\n","      <td>10.213282</td>\n","      <td>1.211617</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["          rank         loo      p_loo     d_loo        weight         se  \\\n","随机斜率模型       0 -171.263208  12.461399  0.000000  6.657864e-01  10.313712   \n","仅包含截距的模型     1 -171.747158  10.166437  0.483949  3.342136e-01  10.156738   \n","随机截距模型       2 -171.944453  10.723415  0.681245  4.440892e-16  10.213282   \n","\n","               dse  warning loo_scale  \n","随机斜率模型    0.000000    False       log  \n","仅包含截距的模型  1.729660    False       log  \n","随机截距模型    1.211617    False       log  "]},"execution_count":31,"metadata":{},"output_type":"execute_result"}],"source":["# 将三个模型的采样结果进行比较\n","compare_dict = {\n","    \"仅包含截距的模型\": trace1, \n","    \"变化截距模型\": trace2, \n","    \"变化斜率模型\": trace3\n","    }\n","# 选择loo方法进行比较\n","comp = az.compare(compare_dict, ic='loo')\n","comp"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"35805896327744FFA347A65339DD54E2","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["结果显示， 变化截距+斜率模型的拟合度好于其他两个模型，这表示，在不同学校间学前教育对留级率的影响是不同的。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"24F66683D1B940EC9DA1AE80C6731AB1","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (9)统计推断"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D5039C8C32C444C3B8AA751FF2B5AB26","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["通过模型比较可以发现，学校所带来的变化效应的影响。如果体现这种影响呐？\n","\n","通过 arviz提供的函数 `plot_forest()` 可以可视化在不同学校间学前教育所带来的影响。\n","- 图中左边的编号是 学校ID。\n","- 94% HDI 是学前教育效应 $\\beta$ 在不同学校间的可信区间。\n","- 通过 94% HDI判断效应是否可信发现，学前教育效应在一些学校间 (比如10418, 20204)存在差异，而在一些学校 (比如10109， 20309)不存在差异。"]},{"cell_type":"code","execution_count":32,"metadata":{"collapsed":false,"id":"8CBA7B60EE754A6F89738114919F85B0","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/18429B97893648348C8DD6B5B495AC3E/rmd0u4d1qe.png\">"],"text/plain":["<Figure size 432x424.8 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["sort_ = trace3.posterior[\"beta\"].mean(dim=[\"chain\", \"draw\"])\n","sort_ = sort_.sortby(sort_, ascending=True)\n","\n","az.plot_forest(\n","    trace3,\n","    var_names=[\"beta\"],\n","    combined=True,\n","    coords={\"school\": sort_[\"school\"]}\n",")\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"5DFD65BA3338455AB8B3020AE519E5B3","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["同样，我们可以通过 `az.summary` 函数查看各参数的具体数值。这里只展示了 beta参数的数值，大家可自行查看 alpha参数的相关统计值。"]},{"cell_type":"code","execution_count":33,"metadata":{"collapsed":false,"id":"046ABA8FE3D845FD85705E381172F450","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mean</th>\n","      <th>sd</th>\n","      <th>hdi_3%</th>\n","      <th>hdi_97%</th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>beta[0]</th>\n","      <td>-0.947</td>\n","      <td>0.713</td>\n","      <td>-2.269</td>\n","      <td>0.263</td>\n","      <td>0.048</td>\n","      <td>0.034</td>\n","      <td>199.0</td>\n","      <td>932.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[1]</th>\n","      <td>-0.176</td>\n","      <td>0.577</td>\n","      <td>-1.257</td>\n","      <td>0.968</td>\n","      <td>0.020</td>\n","      <td>0.015</td>\n","      <td>763.0</td>\n","      <td>1690.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[2]</th>\n","      <td>-0.618</td>\n","      <td>0.558</td>\n","      <td>-1.635</td>\n","      <td>0.535</td>\n","      <td>0.019</td>\n","      <td>0.013</td>\n","      <td>842.0</td>\n","      <td>1817.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>beta[3]</th>\n","      <td>0.085</td>\n","      <td>0.599</td>\n","      <td>-1.102</td>\n","      <td>1.170</td>\n","      <td>0.028</td>\n","      <td>0.020</td>\n","      <td>458.0</td>\n","      <td>1097.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[4]</th>\n","      <td>-1.290</td>\n","      <td>0.835</td>\n","      <td>-2.753</td>\n","      <td>0.138</td>\n","      <td>0.080</td>\n","      <td>0.057</td>\n","      <td>124.0</td>\n","      <td>1198.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[5]</th>\n","      <td>-0.680</td>\n","      <td>0.906</td>\n","      <td>-2.262</td>\n","      <td>1.018</td>\n","      <td>0.094</td>\n","      <td>0.067</td>\n","      <td>100.0</td>\n","      <td>271.0</td>\n","      <td>1.03</td>\n","    </tr>\n","    <tr>\n","      <th>beta[6]</th>\n","      <td>-0.846</td>\n","      <td>0.667</td>\n","      <td>-2.154</td>\n","      <td>0.273</td>\n","      <td>0.040</td>\n","      <td>0.028</td>\n","      <td>284.0</td>\n","      <td>1002.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[7]</th>\n","      <td>-1.079</td>\n","      <td>0.694</td>\n","      <td>-2.310</td>\n","      <td>0.311</td>\n","      <td>0.034</td>\n","      <td>0.024</td>\n","      <td>393.0</td>\n","      <td>1226.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[8]</th>\n","      <td>-0.386</td>\n","      <td>0.701</td>\n","      <td>-1.697</td>\n","      <td>0.982</td>\n","      <td>0.034</td>\n","      <td>0.024</td>\n","      <td>336.0</td>\n","      <td>879.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[9]</th>\n","      <td>-0.888</td>\n","      <td>0.696</td>\n","      <td>-2.199</td>\n","      <td>0.370</td>\n","      <td>0.053</td>\n","      <td>0.038</td>\n","      <td>178.0</td>\n","      <td>1251.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[10]</th>\n","      <td>-0.932</td>\n","      <td>0.660</td>\n","      <td>-2.234</td>\n","      <td>0.251</td>\n","      <td>0.028</td>\n","      <td>0.020</td>\n","      <td>570.0</td>\n","      <td>1433.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[11]</th>\n","      <td>-0.496</td>\n","      <td>0.646</td>\n","      <td>-1.676</td>\n","      <td>0.690</td>\n","      <td>0.053</td>\n","      <td>0.038</td>\n","      <td>138.0</td>\n","      <td>1085.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[12]</th>\n","      <td>-0.723</td>\n","      <td>0.589</td>\n","      <td>-1.931</td>\n","      <td>0.366</td>\n","      <td>0.022</td>\n","      <td>0.016</td>\n","      <td>686.0</td>\n","      <td>1179.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[13]</th>\n","      <td>-0.145</td>\n","      <td>0.521</td>\n","      <td>-1.190</td>\n","      <td>0.798</td>\n","      <td>0.023</td>\n","      <td>0.016</td>\n","      <td>487.0</td>\n","      <td>1561.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[14]</th>\n","      <td>-0.010</td>\n","      <td>0.590</td>\n","      <td>-1.118</td>\n","      <td>1.038</td>\n","      <td>0.047</td>\n","      <td>0.033</td>\n","      <td>160.0</td>\n","      <td>1128.0</td>\n","      <td>1.02</td>\n","    </tr>\n","    <tr>\n","      <th>beta[15]</th>\n","      <td>-0.326</td>\n","      <td>0.740</td>\n","      <td>-1.721</td>\n","      <td>1.153</td>\n","      <td>0.021</td>\n","      <td>0.016</td>\n","      <td>1121.0</td>\n","      <td>1754.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[16]</th>\n","      <td>-0.362</td>\n","      <td>0.552</td>\n","      <td>-1.434</td>\n","      <td>0.673</td>\n","      <td>0.018</td>\n","      <td>0.013</td>\n","      <td>937.0</td>\n","      <td>1664.0</td>\n","      <td>1.01</td>\n","    </tr>\n","    <tr>\n","      <th>beta[17]</th>\n","      <td>-0.573</td>\n","      <td>0.631</td>\n","      <td>-1.741</td>\n","      <td>0.618</td>\n","      <td>0.036</td>\n","      <td>0.026</td>\n","      <td>297.0</td>\n","      <td>1435.0</td>\n","      <td>1.00</td>\n","    </tr>\n","    <tr>\n","      <th>beta[18]</th>\n","      <td>-0.679</td>\n","      <td>0.594</td>\n","      <td>-1.760</td>\n","      <td>0.412</td>\n","      <td>0.041</td>\n","      <td>0.029</td>\n","      <td>204.0</td>\n","      <td>1436.0</td>\n","      <td>1.01</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["           mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n","beta[0]  -0.947  0.713  -2.269    0.263      0.048    0.034     199.0   \n","beta[1]  -0.176  0.577  -1.257    0.968      0.020    0.015     763.0   \n","beta[2]  -0.618  0.558  -1.635    0.535      0.019    0.013     842.0   \n","beta[3]   0.085  0.599  -1.102    1.170      0.028    0.020     458.0   \n","beta[4]  -1.290  0.835  -2.753    0.138      0.080    0.057     124.0   \n","beta[5]  -0.680  0.906  -2.262    1.018      0.094    0.067     100.0   \n","beta[6]  -0.846  0.667  -2.154    0.273      0.040    0.028     284.0   \n","beta[7]  -1.079  0.694  -2.310    0.311      0.034    0.024     393.0   \n","beta[8]  -0.386  0.701  -1.697    0.982      0.034    0.024     336.0   \n","beta[9]  -0.888  0.696  -2.199    0.370      0.053    0.038     178.0   \n","beta[10] -0.932  0.660  -2.234    0.251      0.028    0.020     570.0   \n","beta[11] -0.496  0.646  -1.676    0.690      0.053    0.038     138.0   \n","beta[12] -0.723  0.589  -1.931    0.366      0.022    0.016     686.0   \n","beta[13] -0.145  0.521  -1.190    0.798      0.023    0.016     487.0   \n","beta[14] -0.010  0.590  -1.118    1.038      0.047    0.033     160.0   \n","beta[15] -0.326  0.740  -1.721    1.153      0.021    0.016    1121.0   \n","beta[16] -0.362  0.552  -1.434    0.673      0.018    0.013     937.0   \n","beta[17] -0.573  0.631  -1.741    0.618      0.036    0.026     297.0   \n","beta[18] -0.679  0.594  -1.760    0.412      0.041    0.029     204.0   \n","\n","          ess_tail  r_hat  \n","beta[0]      932.0   1.01  \n","beta[1]     1690.0   1.01  \n","beta[2]     1817.0   1.02  \n","beta[3]     1097.0   1.01  \n","beta[4]     1198.0   1.01  \n","beta[5]      271.0   1.03  \n","beta[6]     1002.0   1.01  \n","beta[7]     1226.0   1.01  \n","beta[8]      879.0   1.01  \n","beta[9]     1251.0   1.01  \n","beta[10]    1433.0   1.01  \n","beta[11]    1085.0   1.01  \n","beta[12]    1179.0   1.01  \n","beta[13]    1561.0   1.01  \n","beta[14]    1128.0   1.02  \n","beta[15]    1754.0   1.01  \n","beta[16]    1664.0   1.01  \n","beta[17]    1435.0   1.00  \n","beta[18]    1436.0   1.01  "]},"execution_count":33,"metadata":{},"output_type":"execute_result"}],"source":["az.summary(trace3, var_names = [\"beta\"])"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"BF2461AF272C417DB33FD4FF2E8A8ED5","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["前面的分析是针对 population level 参数，即变化效应。而我们一般想了解组层面 (group level)的效应是否存在，即固定效应。\n","\n","组层面的效应对应的参数是 mu_alpha, mu_beta, sigma_alpha, sigma_beta。我们可以通过 `az.summary` 与 `az.plot_posterior` 来查看这些参数。"]},{"cell_type":"code","execution_count":36,"metadata":{"collapsed":false,"id":"56387A83EC9649E792DCBC2966C38F4E","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/958BB12B156C4B5882157AF0D23931A1/rmd0xfn072.png\">"],"text/plain":["<Figure size 993.6x331.2 with 2 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_posterior(trace3, var_names=['mu_alpha','mu_beta'])\n","plt.show()"]},{"cell_type":"code","execution_count":37,"metadata":{"collapsed":false,"id":"0A8CCC33CE2847288A9744B7F06C11EB","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mean</th>\n","      <th>sd</th>\n","      <th>hdi_3%</th>\n","      <th>hdi_97%</th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>mu_alpha</th>\n","      <td>-1.014</td>\n","      <td>0.303</td>\n","      <td>-1.613</td>\n","      <td>-0.454</td>\n","      <td>0.014</td>\n","      <td>0.010</td>\n","      <td>429.0</td>\n","      <td>871.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>mu_beta</th>\n","      <td>-0.583</td>\n","      <td>0.394</td>\n","      <td>-1.290</td>\n","      <td>0.196</td>\n","      <td>0.021</td>\n","      <td>0.015</td>\n","      <td>342.0</td>\n","      <td>899.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["           mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n","mu_alpha -1.014  0.303  -1.613   -0.454      0.014    0.010     429.0   \n","mu_beta  -0.583  0.394  -1.290    0.196      0.021    0.015     342.0   \n","\n","          ess_tail  r_hat  \n","mu_alpha     871.0    1.0  \n","mu_beta      899.0    1.0  "]},"execution_count":37,"metadata":{},"output_type":"execute_result"}],"source":["az.summary(trace3, var_names=['mu_alpha','mu_beta'])"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"8F8D5A2C354E4FD78F4BF55ED2B651E2","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["whoops！ 别忘了广义线性模型中链接函数的存在，我们需要把上述的参数值转换到0-1的范围。\n","\n","可以通过 inv.logit函数进行转换，即 $1 / (1 + exp(-\\theta))$。"]},{"cell_type":"code","execution_count":38,"metadata":{"collapsed":false,"id":"74A36E8E283B47FBA413AB31F6D18998","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stdout","output_type":"stream","text":["未接受学前教育学生的留级率(%) =  0.26611387752965954 \n"," 接受学前教育学生的留级率(%) =  0.6244331002322159 \n"," 学前教育对留级率的影响(%) =  0.3583192227025564\n"]}],"source":["# 模型预测值\n","p_PPED0 = 1 / (1 + np.exp(-trace3.posterior[\"mu_alpha\"].mean())).to_pandas()\n","p_PPED1 = 1 / (1 + np.exp(-(trace3.posterior[\"mu_alpha\"].mean()+trace3.posterior[\"mu_beta\"].mean()))).to_pandas()\n","print(\n","    \"未接受学前教育学生的留级率(%) = \",p_PPED0, \n","    \"\\n 接受学前教育学生的留级率(%) = \", p_PPED1,\n","    \"\\n 学前教育对留级率的影响(%) = \", p_PPED0 - p_PPED1\n","    )"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["# 真实数据值\n","p_PPED = data.groupby('PPED').REPEAT.mean() \n","p_PPED0 = p_PPED[0]\n","p_PPED1 = p_PPED[1]\n","print(\n","    \"未接受学前教育学生的留级率(%) = \",p_PPED0, \n","    \"\\n 接受学前教育学生的留级率(%) = \", p_PPED1,\n","    \"\\n 学前教育对留级率的影响(%) = \", p_PPED0 - p_PPED1\n","    )"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D94FD6E6D4504B218394CB671D98A6FA","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["最后，我们比较，模型3(变化截距+斜率)和模型2(变化截距模型)中学前教育组层面的效应是否相同。\n","\n","可见，两个模型的估计非常类似。但是变化截距+斜率模型比变化截距模型更多地考虑了学前教育效应在不同学校间的差异。"]},{"cell_type":"code","execution_count":39,"metadata":{"collapsed":false,"id":"C269B6C1C1F64BC39675B850BFD60162","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stdout","output_type":"stream","text":["随机截距模型中学前教育的效应(%) =  0.3881406175872315 \n"," 随机斜率模型中学前教育的效应(%) =  0.3583192227025564\n"]}],"source":["beta1 = 1 / (1 + np.exp(-(trace2.posterior[\"beta\"].mean()))).to_pandas()\n","beta2 = 1 / (1 + np.exp(-(trace3.posterior[\"mu_beta\"].mean()))).to_pandas()\n","print(\n","    \"变化截距模型中学前教育的效应(%) = \", beta1, \n","    \"\\n 变化截距+斜率模型中学前教育的效应(%) = \", beta2\n","    )"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"B021F312A2A347BBB1207B5728E4A877","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 练习\n","\n","前面我们定义了各种分层模型，但我们没有演示只包含变化斜率的模型，如下图。\n","\n","接下来的练习大家可以尝试以下定义只包含变化斜率的贝叶斯广义线性分层模型。\n","\n","关键在于：\n","- 只允许斜率 beta 随学习变化。而截距 alpha 不变化。\n","- 如何定义分层模型的先验以及相应的超先验。\n","\n","![Image Name](https://cdn.kesci.com/upload/image/rmank74ytm.png?imageView2/0/w/640/h/640)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D4294789781C40588E0D97BFA44332F0","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["结合上节课学习的GLM和HLM，若我们的数据因变量为离散数据，且数据为嵌套结构时，我们可以将HLM和GLM结合。\n","\n"]},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"0E80277E989A412A87374BB073AE0B7B","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["#加载需要使用的库\n","%matplotlib inline\n","import numpy as np \n","from scipy import stats\n","import matplotlib.pyplot as plt\n","import pandas as pd\n","import arviz as az\n","import pymc3 as pm\n","\n","np.random.seed(123)  # 变化数种子，确保随后生成的变化数相同\n","\n","# 加载数据\n","data = pd.read_csv(\"/home/mw/input/data1317/thaieduc.csv\")\n","data.PPED = data.PPED.map({'yes':1,'no':0})\n","\n","# 将数据分层变量\"学校(school)\"转换为因子(factor)类型\n","school_idxs, school = pd.factorize(data.SCHOOLID)\n","# 定义学校与数据的映射：即标注哪名学生(行)属于哪一所学校\n","coords = {\n","    \"school\": school,\n","    \"obs_id\": np.arange(len(school_idxs)),\n","}"]},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"DC1B6FA9EB62440BAA7ECBDA43720232","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model(coords=coords) as model_practice:\n","    ############################################################################\n","    # 要求：定义level2学校的 Hyperpriors\n","    # 提示：通过 mu_beta = pm.Normal(???) 进行定义\n","    ############################################################################\n","    mu_beta = pm.Normal(\"mu_beta\", ...) \n","    sigma_beta = pm.HalfNormal(\"sigma_beta\", ...)  # 注意这是 HalfNormal 分布\n","    \n","    ############################################################################\n","    # 要求：定义先验\n","    # 提示：通过 beta = pm.Normal(???) 进行定义。通过 dims=\"school\" 使得 beta随学校变化\n","    ############################################################################\n","    beta = pm.Normal('beta', ...)\n","    alpha = pm.Normal('alpha', ...)\n","\n","    # 定义数据分层变量\"学校(school)\"\n","    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n","    x = pm.Data(\"x\", data.PPED, dims=\"obs_id\")\n","    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n","    p = pm.Deterministic(\"p\", alpha + beta[school_idx]*x)\n","    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")"]},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"474F5F5FA2634543A708A01FE3C168EE","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["pm.model_to_graphviz(model_practice)"]},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"2ECCFD958A0D487C99DD6A1DA1D78BB9","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with model_practice:\n","    trace_practice = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"F21BB91B5F5C4007BC768C8488F24721","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["az.summary(trace_practice, var_names=['mu_beta'])"]},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"3C98890BDD304E4A820F7A54E48466C8","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["az.plot_posterior(trace_practice, var_names=['mu_beta'])\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"FDCD51C84169478288BBA988269E57F7","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["## Part 3: 示例2——信号检测论"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A8DB0337AF1B4B389B909DE9B1B345E4","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### Workflow\n","\n","![Image Name](https://cdn.kesci.com/upload/image/rkvikqg9q6.png?imageView2/0/w/650/h/650)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C5EADC95006343F78B7CC7A6D534E426","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (1) 提出研究问题"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"FE53C2207CA34D1BBBF594FAAF597992","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["假如我们在进行一项再认实验，这个实验向被试呈现了一组图片，这些图片中有些是新的（被试之前没有见过的），而有些是旧的（被试之前已看过的）。\n","\n","被试需要对呈现的每一张图片进行反应：“新（未见过）”或者“旧（见过）”。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6616AE8C2B1D4A95BD6DAD40F9E20898","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["信号检测论(Signal Detection Theory, SDT)通过分别测量被试的敏感性以及判断标准来建立每个被试再认过程的模型。\n","\n","信号检测论模型的基本理念是：在每个试次中，当刺激呈现时，被试心里会产生一种“熟悉”（或者记忆强度）的信号。随后被试会依据这些信号，判断当前的刺激是“新”或“旧”。\n","![Image Name](https://th.bing.com/th/id/R.1565c09f95a7411070316c78242814c8?rik=rKBAJumpMrZ5Qg&riu=http%3a%2f%2fmedia.open.com.cn%2fmedia_file004%2f1003%2fdongshi%2fshiyanxlx%2fchapter6%2fimage%2f01.jpg&ehk=ui4XLQbTaXgMxP3BuNoBRAnww2cCrMLys2tN%2f1vm%2fvA%3d&risl=&pid=ImgRaw&r=0)\n","\n","被试的反应有信号检测论的参数d'（辨别力指标/感觉敏感性指标）和 c（报告标准）决定。\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"194E35BF46FF4638B366D710E69C4F36","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["那么我们该如何计算被试的参数d'（辨别力指标/感觉敏感性指标）和 c（报告标准）？"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"73A15DDC901B4F00AF9B607DF93DC929","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (2) 数据收集"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"25D4EA190C044351899224E2EABA66C3","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["我们使用Skagerberg和Wright记忆一致性(memory conformity)研究中的对照组(Skagerberg & Wright，2008)"]},{"cell_type":"code","execution_count":3,"metadata":{"collapsed":false,"id":"D7E7E3B7982F460F9F9F33F669070729","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>Unnamed: 0</th>\n","      <th>subno</th>\n","      <th>sayold</th>\n","      <th>isold</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>1</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>0</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>2</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>1</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>3</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>1</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>4</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>1</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>5</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["   Unnamed: 0  subno  sayold  isold\n","0           1     53       1      0\n","1           2     53       1      1\n","2           3     53       1      1\n","3           4     53       1      1\n","4           5     53       1      0"]},"execution_count":3,"metadata":{},"output_type":"execute_result"}],"source":["# 加载数据\n","data = pd.read_csv(\"/home/mw/input/data1317/confcontr.csv\")\n","# 查看数据\n","data.head()\n","# 该数据中每一行都是一个试次，subno为被试编号，sayold为被试的新旧判断，isold为图片本身的新旧程度"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"7D5DB2D697A843FAA7EF5A0495496940","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 描述统计与探索性分析"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F39E38E50B334AA3892C02DC5ACD8DD4","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["首先，我们可以查看单个被试每个试次的反应（击中率(hit)、虚报率(false alarm)、正确拒斥率(correct rejection)、漏报率(miss)）"]},{"cell_type":"code","execution_count":55,"metadata":{"collapsed":false,"id":"EBAC45D04E7D4DD68E518391F1866EBE","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stdout","output_type":"stream","text":["hr is 0.5319148936170213 \n","mr is 0.46808510638297873 \n","crr is 0.6226415094339622 \n","far is 0.37735849056603776 \n"]}],"source":["data = data.loc[data['subno']==53]\n","hr = len(data.loc[(data['sayold']==1)&(data['isold']==1)])/len(data.loc[data['isold']==1])\n","mr = len(data.loc[(data['sayold']==0)&(data['isold']==1)])/len(data.loc[data['isold']==1])\n","crr = len(data.loc[(data['sayold']==0)&(data['isold']==0)])/len(data.loc[data['isold']==0])\n","far = len(data.loc[(data['sayold']==1)&(data['isold']==0)])/len(data.loc[data['isold']==0])\n","print(f\"hr is {hr} \")\n","print(f\"mr is {mr} \")\n","print(f\"crr is {crr} \")\n","print(f\"far is {far} \")"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"EBB6CD49FF314626BD71E05C7D14B825","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (3) 选择模型"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9EBAC547DD354452962182530410A4CF","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["最常见的SDT模型，即假定被试对两种信号的熟悉程度的分布属于高斯分布（即正态分布），且这两个分布的方差是齐性的，但均值可能不同（例如，对旧刺激的熟悉程度更高）。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"2D50B22803F7482BA06CC52B5054449D","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["在方差齐性的情况下，对于单个被试，d'可由标准化后的击中率和虚报率(hit and false alarm rates)之差计算得来(Stanislaw & Todorov，1999)\n","\n","$d'=\\phi^{-1}(HR)-\\phi^{-1}(FAR)$\n","\n","$\\phi$是累积正态密度函数(cumulative normal density function)，可将z分数转换为概率，其反函数$\\phi^{-1}$ 将概率（例如击中率或虚报率）转换为z分数。标准化后的击中率和虚报率分别表示为 $z(HR)$ 和$z(FAR)$。\n","\n","反应标准$c=-z(FAR)$（负虚报率的标准分数）表示(DeCarlo，1998)。\n","\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"03F8A48F5C1241BEA3EAD758FCFBCF04","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["在除了直接对d'和c进行计算，我们可以使用广义线性模型来进行计算。\n","\n","在二项反应的GLM中，我们假定反应属于Bernoulli分布.\n","\n","即：$y_i \\sim Bernoulli(p_i)$"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F0478F0E9ED84A3497FB520CF32AA300","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["被试反应即对图片新旧程度的判断是由击中率(HR)或虚报率(FAR)来决定的，故概率p表示击中率(HR)。\n","\n","概率p本身的范围为0-1，因此，我们不能直接用线性模型直接对p进行建模。\n","\n","相反，我们用一个连结函数将p转化为‘线性预测值’ $η$，并用线性回归模型对$η$建模。\n","\n","即：$p_i = \\phi(\\eta_i)$\n","\n","如前所述，$\\phi$是一个累积正态密度函数并能将z 值转化为概率，那我们可以用一个截距和斜率来表示$\\eta$。\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"0E753B3E088D4C458538F0D5BEC0C025","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["之所以用$\\phi$而不是logit函数，是因为击中率或虚报率本身是一个概率，它是由信号分布的z分数转化而来，而不是简单地将至于范围转化到[0,1]\n","\n","即：$\\eta_i = \\beta_0 + \\beta_1 * isold_i$\n","\n","由于$d'=\\phi^{-1}(HR)-\\phi^{-1}(FAR) = z(HR) - z(FAR)$，$c=-z(FAR)$\n","\n","当isold为0，反应由虚报率决定，即$z(FAR) = -c$\n","\n","当isold为1，反应由击中率决定，即$z(HR) = z(HR) - z(FAR) + z(FAR) = -c+d$\n","\n","故在线性模型中，$\\beta_0$代表c(报告标准)的相反数，$\\beta_1$代表d(辨别力指标)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"288BFB9F301E49AC852AAB594EA1F31F","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["我们首先展示一个单被试模型的例子，由于该模型只有一个被试，故是一个单层的模型\n","$$\n","y_{(i)} \\sim Bernoulli(p_{(i)}) \\\\\n","p_{(i)} = \\phi(\\eta) \\\\\n","\\eta = -c + d *x_{isold}\n","$$\n","\n","\n","**先验**\n","\n","$$\n","\\ c \\sim Normal(0,1) \\\\\n","\\ d \\sim Normal(0,1) \\\\ \n","$$"]},{"cell_type":"code","execution_count":59,"metadata":{"collapsed":false,"id":"B26C0961744A43B5AEFADD4E5B6108BF","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model(coords=coords) as model0:\n","    # 定义priors\n","    c = pm.Normal(\"c\", mu=0.0, sigma=1)  \n","    d = pm.Normal(\"d\", mu=0.0, sigma=1)\n","    # 自变量x是之前已经载入的数据\n","    x = pm.Data(\"x\", data['isold'])\n","    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n","    eta = pm.Deterministic(\"eta\", -c + d*x)\n","    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    sayold = pm.Bernoulli(\"syaold\", logit_p=eta, observed=data['sayold'])"]},{"cell_type":"code","execution_count":60,"metadata":{"collapsed":false,"id":"1EA98DD253F44485A5BA37853E6A346A","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/DC9DF293D73E4296A4710638CF82F673/rmd94k89uk.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7fe1f46e0550>"]},"execution_count":60,"metadata":{},"output_type":"execute_result"}],"source":["pm.model_to_graphviz(model0)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"735E2DE057264A5FBC6092A7AF924B9B","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["单被试模型仅为简单说明如何用GLM建立信号检测论，接下来的过程为了避免混淆将不对该模型做进一步的分析"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"723D322B02564238961CF360E0CE71A3","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["接下来我们将展示在分层结构中，我们如何利用GLMM建立信号检测论的模型"]},{"cell_type":"code","execution_count":74,"metadata":{"collapsed":false,"id":"8B3525194B0046B8909682ACD31AA9A2","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>Unnamed: 0</th>\n","      <th>subno</th>\n","      <th>sayold</th>\n","      <th>isold</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>1</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>0</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>2</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>1</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>3</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>1</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>4</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>1</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>5</td>\n","      <td>53</td>\n","      <td>1</td>\n","      <td>0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["   Unnamed: 0  subno  sayold  isold\n","0           1     53       1      0\n","1           2     53       1      1\n","2           3     53       1      1\n","3           4     53       1      1\n","4           5     53       1      0"]},"execution_count":74,"metadata":{},"output_type":"execute_result"}],"source":["data = pd.read_csv(\"/home/mw/input/data1317/confcontr.csv\")\n","data.head()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"5912E0CB26564DA9AB440A552274AC15","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 模型1：变化截距模型\n","\n","我们假定各被试的c(报告标准)是不同的\n","\n","**固定效应**\n","$$\n","y_{i,j} \\sim Bernoulli(p_{i,j}) \\\\\n","p_{i,j} = \\phi(\\eta_{i,j}) \\\\\n","\\eta_{i,j} = -c_{i,j} + d_{i,j} *x_{isold,i,j}\n","$$\n","\n","\n","**先验**\n","\n","$$\n","\\ c_{i,j} \\sim Normal(\\mu_{c,j},\\sigma_{c,j}) \\\\\n","\\ d_{i,j} \\sim Normal(0,1) \\\\ \n","$$\n","**随机效应**\n","$$\n","\\ \\mu_{c,j} \\sim Normal(0,1) \\\\\n","\\ \\sigma_{c,j} \\sim HalfNormal(1) \\\\\n","$$"]},{"cell_type":"code","execution_count":87,"metadata":{"collapsed":false,"id":"368D5EF3C5174D3A8217382A3BC4186F","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["# 将数据分层变量\"被试(subject)\"转换为因子(factor)类型\n","subj_idxs, subject = pd.factorize(data.subno)\n","# 定义学校与数据的映射：即标注哪名学生(行)属于哪一所学校\n","coords = {\n","    \"subject\": subject,\n","    \"obs_id\": np.arange(len(subj_idxs)),\n","}"]},{"cell_type":"code","execution_count":88,"metadata":{"collapsed":false,"id":"4B2092910D5B4CE6A43186AE133C2E66","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model(coords=coords) as model1:\n","    # 定义level2被试的 Hyperpriorspriors\n","    mu_c = pm.Normal(\"mu_c\", mu=0.0, sigma=1)  \n","    sigma_c = pm.HalfNormal(\"sigma_c\", sigma = 1)\n","    # 定义priors\n","    c = pm.Normal(\"c\", mu=mu_c, sigma=sigma_c,dims=\"subject\")  \n","    d = pm.Normal(\"d\", mu=0.0, sigma=1)\n","    # 定义数据分层变量\"被试(subject)\"\n","    subj_idx = pm.Data(\"subject\", subj_idxs, dims=\"obs_id\")\n","    # 自变量x是之前已经载入的数据\n","    x = pm.Data(\"x\", data['isold'], dims=\"obs_id\")\n","    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n","    eta = pm.Deterministic(\"eta\", -c[subj_idx] + d*x)\n","    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    sayold = pm.Bernoulli(\"syaold\", logit_p=eta, observed=data['sayold'],dims=\"obs_id\")"]},{"cell_type":"code","execution_count":89,"metadata":{"collapsed":false,"id":"4B4A02B41EBD435EA3A48AC106BD2CC0","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/E99851964F454E04B49689C61C0FFB0B/rmdarygr9c.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7fe164e4b510>"]},"execution_count":89,"metadata":{},"output_type":"execute_result"}],"source":["pm.model_to_graphviz(model1)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"2401BD6434904450A19A99A7ED38E152","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 模型2：变化斜率模型\n","\n","我们假定各被试的d(辨别力指标)是不同的\n","\n","**固定效应**\n","$$\n","y_{i,j} \\sim Bernoulli(p_{i,j}) \\\\\n","p_{i,j} = \\phi(\\eta_{i,j}) \\\\\n","\\eta_{i,j} = -c_{i,j} + d_{i,j} *x_{isold,i,j}\n","$$\n","\n","\n","**先验**\n","\n","$$\n","\\ c_{i,j} \\sim Normal(0,1) \\\\\n","\\ d_{i,j} \\sim Normal(\\mu_{d,j},\\sigma_{d,j}) \\\\ \n","$$\n","**随机效应**\n","$$\n","\\ \\mu_{d,j} \\sim Normal(0,1) \\\\\n","\\ \\sigma_{d,j} \\sim HalfNormal(1) \\\\\n","$$"]},{"cell_type":"code","execution_count":92,"metadata":{"collapsed":false,"id":"D4C2CD117DA24E4A9988AC2EF9607E13","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model(coords=coords) as model2:\n","    # 定义level2被试的 Hyperpriorspriors\n","    mu_d = pm.Normal(\"mu_d\", mu=0.0, sigma=1)  \n","    sigma_d = pm.HalfNormal(\"sigma_d\", sigma = 1)\n","    # 定义priors\n","    c = pm.Normal(\"c\", mu=0, sigma=1)  \n","    d = pm.Normal(\"d\", mu=mu_d, sigma=sigma_d, dims=\"subject\")\n","    # 定义数据分层变量\"被试(subject)\"\n","    subj_idx = pm.Data(\"subject\", subj_idxs, dims=\"obs_id\")\n","    # 自变量x是之前已经载入的数据\n","    x = pm.Data(\"x\", data['isold'], dims=\"obs_id\")\n","    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n","    eta = pm.Deterministic(\"eta\", -c + d[subj_idx]*x)\n","    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    sayold = pm.Bernoulli(\"syaold\", logit_p=eta, observed=data['sayold'],dims=\"obs_id\")"]},{"cell_type":"code","execution_count":93,"metadata":{"collapsed":false,"id":"488EA150856D49D6B89AEA0B251FA861","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/A4F0B9D023EC4A33A2FAB528BC49EC02/rmdb3c7e7t.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7fe1f4db0f10>"]},"execution_count":93,"metadata":{},"output_type":"execute_result"}],"source":["pm.model_to_graphviz(model2)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"174EBD43169B4B1BA9BCB7302E608A55","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (4)选择先验"]},{"cell_type":"code","execution_count":103,"metadata":{"collapsed":false,"id":"7C81C72F627C45E78A6F6632B2C2CF46","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with model1:\n","    # 先验预测检查\n","    prior_checks = pm.sample_prior_predictive(samples=500)"]},{"cell_type":"code","execution_count":104,"metadata":{"collapsed":false,"id":"24A1BDFAED144DB2A729F2F2EBDEDE87","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/5A6015A47ADB44A88001546D2782FB91/rmdbd4e6q0.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_density(\n","    {'eta':prior_checks['eta']}\n","    )\n","plt.show()"]},{"cell_type":"code","execution_count":108,"metadata":{"collapsed":false,"id":"FB7BF705491E4B68AB47A9694AD0F17F","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with model2:\n","    # 先验预测检查\n","    prior_checks = pm.sample_prior_predictive(samples=500)"]},{"cell_type":"code","execution_count":109,"metadata":{"collapsed":false,"id":"44012BAA026B4698864CA31AE21A1B06","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/A09076FA9A574C6583D2B299A40466E5/rmdbee6o6c.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_density(\n","    {'eta':prior_checks['eta']}\n","    )\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F3CFEEE172A14D1B9D1059709DF75C10","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (5) 拟合数据"]},{"cell_type":"code","execution_count":110,"metadata":{"collapsed":false,"id":"9D42BA673A064DD9BF0517055CF6120F","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [d, c, sigma_c, mu_c]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 17 seconds.\n"]},{"name":"stdout","output_type":"stream","text":["0, dim: obs_id, 3100 =? 3100\n"]}],"source":["with model1:\n","    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n","    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n","    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n","    trace1 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"code","execution_count":111,"metadata":{"collapsed":false,"id":"2BB76A7C08984F9A8CB116A88365F8D0","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [d, c, sigma_d, mu_d]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 16 seconds.\n"]},{"name":"stdout","output_type":"stream","text":["0, dim: obs_id, 3100 =? 3100\n"]}],"source":["with model2:\n","    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n","    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n","    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n","    trace2 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"5984662B1E154590AC732AB241EC74D2","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (6)采样过程诊断\n","\n","如果使用MCMC对后验进行近似，则需要首先对MCMC过程进行评估。\n","\n","* 是否收敛；\n","* 是否接近真实的后验。\n","\n","对采样过程的评估我们会采用目视检查或rhat这个指标"]},{"cell_type":"code","execution_count":112,"metadata":{"collapsed":false,"id":"142944CE11A64B619842ED7BE78F9912","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mean</th>\n","      <th>sd</th>\n","      <th>hdi_3%</th>\n","      <th>hdi_97%</th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>mu_c</th>\n","      <td>1.052</td>\n","      <td>0.084</td>\n","      <td>0.894</td>\n","      <td>1.208</td>\n","      <td>0.002</td>\n","      <td>0.001</td>\n","      <td>2121.0</td>\n","      <td>2562.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>sigma_c</th>\n","      <td>0.330</td>\n","      <td>0.067</td>\n","      <td>0.213</td>\n","      <td>0.459</td>\n","      <td>0.002</td>\n","      <td>0.001</td>\n","      <td>1690.0</td>\n","      <td>1575.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>d</th>\n","      <td>1.677</td>\n","      <td>0.080</td>\n","      <td>1.521</td>\n","      <td>1.825</td>\n","      <td>0.002</td>\n","      <td>0.001</td>\n","      <td>1805.0</td>\n","      <td>2271.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["          mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n","mu_c     1.052  0.084   0.894    1.208      0.002    0.001    2121.0   \n","sigma_c  0.330  0.067   0.213    0.459      0.002    0.001    1690.0   \n","d        1.677  0.080   1.521    1.825      0.002    0.001    1805.0   \n","\n","         ess_tail  r_hat  \n","mu_c       2562.0    1.0  \n","sigma_c    1575.0    1.0  \n","d          2271.0    1.0  "]},"execution_count":112,"metadata":{},"output_type":"execute_result"}],"source":["az.summary(trace1, var_names = [\"mu_c\",\"sigma_c\",\"d\"])"]},{"cell_type":"code","execution_count":115,"metadata":{"collapsed":false,"id":"8C6D49562B124D2FBC2739D99A7E0297","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mean</th>\n","      <th>sd</th>\n","      <th>hdi_3%</th>\n","      <th>hdi_97%</th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>mu_d</th>\n","      <td>1.659</td>\n","      <td>0.119</td>\n","      <td>1.430</td>\n","      <td>1.873</td>\n","      <td>0.003</td>\n","      <td>0.002</td>\n","      <td>1678.0</td>\n","      <td>2299.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>sigma_d</th>\n","      <td>0.492</td>\n","      <td>0.099</td>\n","      <td>0.320</td>\n","      <td>0.688</td>\n","      <td>0.002</td>\n","      <td>0.002</td>\n","      <td>1592.0</td>\n","      <td>1679.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>c</th>\n","      <td>1.032</td>\n","      <td>0.055</td>\n","      <td>0.922</td>\n","      <td>1.131</td>\n","      <td>0.001</td>\n","      <td>0.001</td>\n","      <td>1754.0</td>\n","      <td>2796.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["          mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n","mu_d     1.659  0.119   1.430    1.873      0.003    0.002    1678.0   \n","sigma_d  0.492  0.099   0.320    0.688      0.002    0.002    1592.0   \n","c        1.032  0.055   0.922    1.131      0.001    0.001    1754.0   \n","\n","         ess_tail  r_hat  \n","mu_d       2299.0    1.0  \n","sigma_d    1679.0    1.0  \n","c          2796.0    1.0  "]},"execution_count":115,"metadata":{},"output_type":"execute_result"}],"source":["az.summary(trace2, var_names = [\"mu_d\",\"sigma_d\",\"c\"])"]},{"cell_type":"code","execution_count":116,"metadata":{"collapsed":false,"id":"AA0BF50732E34F2DABE92C892D98BA01","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/D2C62244A87E413A8481C2CE9721057B/rmdbimwgj7.png\">"],"text/plain":["<Figure size 864x432 with 6 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_trace(trace1, var_names = [\"mu_c\",\"sigma_c\",\"d\"])\n","plt.show()"]},{"cell_type":"code","execution_count":117,"metadata":{"collapsed":false,"id":"84D044BC4DFA460B8FE8BC020C1CED94","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/B211E4BFF7BE4ED5A1B8837C9423BB74/rmdbiyajw2.png\">"],"text/plain":["<Figure size 864x432 with 6 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_trace(trace2, var_names = [\"mu_d\",\"sigma_d\",\"c\"])\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E138385B0422471D8BA8BC0360FC27A4","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (7)模型诊断\n","\n","在MCMC有效的前提下，需要继续检验模型是否能够较好地拟合数据。\n","\n","我们会使用后验预测分布通过我们得到的参数生成一批模拟数据，并将其与真实数据进行对比。"]},{"cell_type":"code","execution_count":118,"metadata":{"collapsed":false,"id":"5F5AD963BEAA4FCD89B87FE3F98EE067","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["\n","    <div>\n","      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n","      100.00% [4000/4000 00:03&lt;00:00]\n","    </div>\n","    "],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"update_display_data"},{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n","  FutureWarning,\n"]}],"source":["# 后验预测分布的计算仍在容器中进行\n","with model1:\n","    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n","    ppc_y = pm.sample_posterior_predictive(trace1.posterior) \n","#将ppc_y转化为InferenceData对象合并到trace中\n","az.concat(trace1, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)"]},{"cell_type":"code","execution_count":120,"metadata":{"collapsed":false,"id":"2DAAD40B722C43C5ABE1AFFA6FC9D7E5","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["\n","    <div>\n","      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n","      100.00% [4000/4000 00:03&lt;00:00]\n","    </div>\n","    "],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"update_display_data"}],"source":["# 后验预测分布的计算仍在容器中进行\n","with model2:\n","    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n","    ppc_y = pm.sample_posterior_predictive(trace2.posterior) \n","#将ppc_y转化为InferenceData对象合并到trace中\n","az.concat(trace2, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)"]},{"cell_type":"code","execution_count":121,"metadata":{"collapsed":false,"id":"0E77A4DC9B374BDEBA7AC13E2B10EBAA","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n","  fig.canvas.print_figure(bytes_io, **kw)\n"]},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/ABBB3C253E3042838B4DE34EC93C1978/rmdbmwx2g.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制后验预测分布\n","az.plot_ppc(trace1)\n","plt.show()"]},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"FC657B6EE1EF449AAA2B56054B658856","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["# 绘制后验预测分布\n","az.plot_ppc(trace2)\n","plt.show()"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E3BC6712B26F45E390F6B64E54FA8784","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (8)模型比较"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"B5F7907292BC42CEB06D08E7CC1E733C","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["在模型比较阶段，我们可以同时比较2个模型：\n","- 模型1：变化截距模型。\n","- 模型2：变化截距+斜率模型。"]},{"cell_type":"code","execution_count":122,"metadata":{"collapsed":false,"id":"549C6ACA71004AA99545D3695ADBD7B1","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n","  \"The default method used to estimate the weights for each model,\"\n"]},{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>rank</th>\n","      <th>loo</th>\n","      <th>p_loo</th>\n","      <th>d_loo</th>\n","      <th>weight</th>\n","      <th>se</th>\n","      <th>dse</th>\n","      <th>warning</th>\n","      <th>loo_scale</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>随机斜率模型</th>\n","      <td>0</td>\n","      <td>-1868.771091</td>\n","      <td>23.095884</td>\n","      <td>0.00000</td>\n","      <td>0.538192</td>\n","      <td>22.228239</td>\n","      <td>0.000000</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","    <tr>\n","      <th>随机截距模型</th>\n","      <td>1</td>\n","      <td>-1870.007642</td>\n","      <td>22.829657</td>\n","      <td>1.23655</td>\n","      <td>0.461808</td>\n","      <td>22.199479</td>\n","      <td>6.045258</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["        rank          loo      p_loo    d_loo    weight         se       dse  \\\n","随机斜率模型     0 -1868.771091  23.095884  0.00000  0.538192  22.228239  0.000000   \n","随机截距模型     1 -1870.007642  22.829657  1.23655  0.461808  22.199479  6.045258   \n","\n","        warning loo_scale  \n","随机斜率模型    False       log  \n","随机截距模型    False       log  "]},"execution_count":122,"metadata":{},"output_type":"execute_result"}],"source":["# 将2个模型的采样结果进行比较\n","compare_dict = {\n","    \"变化截距模型\": trace1, \n","    \"变化截距+斜率模型\": trace2\n","    }\n","# 选择loo方法进行比较\n","comp = az.compare(compare_dict, ic='loo')\n","comp"]},{"cell_type":"code","execution_count":123,"metadata":{"collapsed":false,"id":"C502BE249FA74068BD60251DB6BDAB24","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n","  \"The default method used to estimate the weights for each model,\"\n"]},{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>rank</th>\n","      <th>waic</th>\n","      <th>p_waic</th>\n","      <th>d_waic</th>\n","      <th>weight</th>\n","      <th>se</th>\n","      <th>dse</th>\n","      <th>warning</th>\n","      <th>waic_scale</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>随机斜率模型</th>\n","      <td>0</td>\n","      <td>-1868.758226</td>\n","      <td>23.083019</td>\n","      <td>0.000000</td>\n","      <td>0.538506</td>\n","      <td>22.227931</td>\n","      <td>0.000000</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","    <tr>\n","      <th>随机截距模型</th>\n","      <td>1</td>\n","      <td>-1870.006190</td>\n","      <td>22.828205</td>\n","      <td>1.247964</td>\n","      <td>0.461494</td>\n","      <td>22.199454</td>\n","      <td>6.044934</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["        rank         waic     p_waic    d_waic    weight         se       dse  \\\n","随机斜率模型     0 -1868.758226  23.083019  0.000000  0.538506  22.227931  0.000000   \n","随机截距模型     1 -1870.006190  22.828205  1.247964  0.461494  22.199454  6.044934   \n","\n","        warning waic_scale  \n","随机斜率模型    False        log  \n","随机截距模型    False        log  "]},"execution_count":123,"metadata":{},"output_type":"execute_result"}],"source":["# 将2个模型的采样结果进行比较\n","compare_dict = {\n","    \"变化截距模型\": trace1, \n","    \"变化截距+斜率模型\": trace2\n","    }\n","# 选择loo方法进行比较\n","comp = az.compare(compare_dict, ic='waic')\n","comp"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"0331CD2A16AB4AEF8F304C20F0660F29","jupyter":{},"notebookId":"638d47efbd63e7d181de5a85","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 总结\n","\n","- 本节课学习了分层线性模型的基本概念：包括变化截距与变化斜率、固定效应与变化效应、group level与population level。\n","- 通过 PyMc 分别拟合变化截距模型与变化斜率模型(包含变化截距)。"]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.5.2"}},"nbformat":4,"nbformat_minor":0}
